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INTRODUCTION 


The  purpose  of  this  project  is  to  develop  a  computer  program  that 
will  design  and  implement  a  recursive  digital  filter  to  meet  user  speci¬ 
fications  given  as  inputs  to  the  program.  The  filter  type  is  selected 
from  three  well-known  analog  filter  approximations:  (1)  Butterworth , 

maximally  flat;  (2)  Chebyshev,  equiripple  in  the  passband;  and  (3) 
elliptic  (or  Cauer)  ,  equiripple  in  both  the  passband  and  the  stop 
band.  For  a  low-pass  filter,  for  example,  the  user  specifies  the  filter 
type,  the  filter  order,  the  cutoff  frequency,  and  the  sampling  rate. 
High-pass,  band-pass,  and  band-stop  filters  also  may  be  chosen  for 
design. 

The  program  calculates  the  poles  and  the  zeros  of  the  low-pass 
filter  in  the  analog  s  domain.  The  bilinear  transformation  technique  is 
then  used  to  determine  the  coefficients  of  the  transfer  function  in  the 
z  domain.  If  a  high-pass,  band-pass,  or  band-stop  filter  is  required, 
then  the  coefficients  are  determined  by  a  transformation  in  the  z  domain 
of  the  prototype  low-pass  digital  filter  whose  coefficients  have  been 
determined  from  an  analog  low-pass  filter  using  the  bilinear  transforma¬ 
tion. 


After  determining  the  coefficients  of  the  digital  filter  transfer 
function,  a  plot  of  the  magnitude  and  the  phase  may  be  obtained  if 
desired  by  the  user.  A  list  of  the  coefficients  also  is  printed. 
Another  option  is  to  implement  the  *■  1  ter  as  a  system  of  linear  differ¬ 
ence  equations  and  to  process  input  data  consisting  of  sample  values  at 
any  given  sample  rate  through  the  filter.  The  filtered  output  data 
points  are  then  printed  as  output. 

In  summary,  the  program  performs  the  following: 

a.  Calculates  and  prints  the  coefficients  of  the  desired  digital 
f i 1  ter . 

b.  Plots  the  magnitude  and  phase  responses. 

c.  Implements  the  digital  filter  and  processes  data  through  the 
filter . 

d.  Outputs  the  filtered  data. 


2 .  APPROACH 

There  are  two  general  types  of  digital  filters,  recursive  and  non¬ 
recursive.  The  output  response  of  a  nonrecursive  filter  is  a  function 


of  only  the  present  and  past  values  of  the  input  excitation.  A  recur¬ 
sive  filter  is  one  in  which  the  present  output  response  is  a  function  of 
the  present  and  past  values  of  the  input ,  as  well  as  past  values  of  the 
output . 

As  in  analog  filters,  the  approximation  step  in  the  design  of  digi¬ 
tal  filters  is  the  process  whereby  a  realizable  transfer  function  satis¬ 
fying  prescribed  specifications  is  obtained.  To  be  realizable  as  a 
recursive  filter,  a  transfer  function  must  satisfy  the  following  con¬ 
straints  . 


a.  It  must  be  a  rational  function  of  z  with  real  coefficients. 

b.  Its  poles  must  lie  within  the  unit  circle  of  the  z  plane. 

c.  The  degree  of  the  numerator  polynomial  must  be  equal  to  or  less 
than  that  of  the  denominator  polynomial. 

Recursive  filter  approximations  can  be  obtained  from  analog  filter 
approximations  by  using  the  following  methods,  all  of  which  satisfy  the 
above  constraints: 

Invariant- impulse-response  method 
Matched  z  transformation 
Bilinear  transformation 

Each  method  has  certain  advantages  and  disadvantages. 

2. 1  Invariant- Impulse-Response  Method 

Given  an  analog  filter  transfer  function,  HA(s),  the  invariant- 
impulse-response  method  is  implemented  as  follows: 

a.  Obtain  the  impulse  response  of  analog  filter  H,(t). 

b.  Replace  t  by  nT  in  Hft(t). 

c.  Form  the  z  transform  of  Hft(nT).  This  gives  HD(z),  the 
digital  filter  transfer  function.  T  is  the  sample  period  =  1/fg. 

This  method  gives  good  results  if  HA(jw)  »  0  for  w  >  wg/2. 
However,  aliasing  errors  tend  to  restrict  this  method  to  the  design  of 
all  pole  filters. 

2.2  Matched  z  Transformation 


Given  continuous-time  transfer  function 


Ho  n  (s  -  si) 

H  (s)  =  _± - 

A  N 

.11  (S  -  Pi) 
1=1 


a  corresponding  digital  transfer  function  can  be  formed  as 


Vz) 


(2  +  1)L 


where  L  is  an  integer  equal  to  the  number  of  zeros  at  s  =  00  in  Hft(s). 
This  method  gives  reasonable  results  for  high-pass  and  band-stop  fil¬ 
ters,  although  it  tends  to  distort  the  passb.md  ripple  in  Chebyshev  and 
elliptic  filters.  For  low-pass  and  band-pass  filters,  better  approxima¬ 
tions  can  be  obtained  by  using  the  modified  invariant-impulse-response 
method. 1 


2.3  Bilinear  Transformation 


The  bilinear  transformation  method  yields  a  digital  filter  with 
approximately  the  same  time-domain  response  as  the  original  analog 
filter  for  any  excitation. 


H  ( z)  =  H  (s) 
D  A 


s 


2  z  -  1 
T  z  +  1 


The  bilinear 

unit  circle 

b. 

c . 


transformation  maps  these : 

The  open  right-half  s  plane  onto  the  region  exterior  to  the 
|z|  =  1  of  the  z  plane 

The  j  axis  of  the  s  plane  onto  the  unit  circle  |z|  =1 
The  open  left-half  s  plane  onto  the  interior  of  the  unit 


circle 


From  property  c  it  follows  that  a  stable  analog  filter  yields  a 
stable  digital  filter  and,  since  the  transformation  has  real  coeffi¬ 
cients,  Hd(z)  has  real  coefficients. 


1A.  Antoniou,  Digital  Filters:  Analysis  and  Design,  McGraw-Hill  Book 
Co.,  Inc.,  New  York  (1979). 
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8,(6^ 


If  only  the  amplitude  response  is  of  concern,  the  warping 
effect  can  for  all  practical  purposes  be  eliminated  by  prewarping  the 
analog  filter.  For  example,  if  a  low-pass  cutoff  frequency  of  is 

desired  for  the  digital  filter,  then  the  analog  filter  is  first  designed 
for  an  unnormalized  cutoff  frequency  given  by 


The  phase  response  of  the  derived  digital  filter  is  nonlinear 
because  of  the  warping  effect.  Furthermore ,  little  can  be  done  to 
linearize  it  except  by  employing  delay  equalization.  Consequently,  if 
it  is  mandatory  to  preserve  a  linear  phase  response,  alternative  methods 
should  be  considered. 


The  bilinear  transformation  is  the  most  important  of  the  tech¬ 
niques  used  to  obtain  digital  recursive  filters  from  analog  filters. 
The  passband  and  loss  character  istics  of  the  analog  filter  are  pre¬ 
served,  and  there  is  no  aliasing  effect.  Frequency  distortion  is  com¬ 
pensated  at  and  the  transfer  function  can  be  obtained  by  a 

relatively  easy  transformation.  For  these  reasons,  the  bilinear  trans¬ 
formation  was  chosen  as  the  method  of  obtaining  the  digital  filter 
transfer  function  for  this  project. 


2.4  Analog  Filter  Designs 


The  theory  of  analog  filter  approximations  has  been  extensively 
developed. ^  4 


2. 4. 1  Butterworth 


The  Butterworth  approximation  is  the  simplest  all-pole  type. 
In  the  stop  band,  the  attenuation  characteristics  monotonically  decrease 
as  a  function  of  frequency.  The  Butterworth  filter  is  termed  also  a 
maximally  flat  magnitude  approximation  in  that  the  error  in  the  passband 
is  also  a  monotonically  decreasing  function.  The  equations  for  calcula¬ 
ting  the  pole  locations  are  given  in  appendix  A. 


Antoniou,  Digital  Filters:  Analysis  and  Design,  McGraw-Hill  Book 
Co.,  Inc.,  New  York  (1979). 

2M.  S.  Ghausi,  Principles  and  Design  of  Linear  Active  Circuits, 
McGraw-Hill  Book  Co.,  Inc.,  New  York  (1965),  ch.  4. 

H .  Y.-F.  Lam,  Analog  and  Digital  Filters,  Prentice-Hall,  Inc.,  Engle¬ 
wood  Cliffs,  NJ  (1979). 

4 D .  E.  Johnson,  Introduction  to  Filter  Theory,  Prentice-Hall ,  Inc., 
Englewood  Cliffs,  NJ  (1976). 
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2.4.2  Chebyshev 

A  second  approximation  that  improves  on  the  rate  of  change  of 
the  attenuation  between  passband  and  stop  band  over  that  of  the  Butter- 
worth  filter  is  the  Chebyshev  filter.  The  error  in  the  passband  is 
distributed  evenly  in  an  oscillating  manner.  This  is  called  an  equi- 
ripple  approximation.  In  the  stop  band,  the  magnitude  decreases  mono- 
tonically  with  a  faster  cutoff  rate  than  that  of  a  Butterworth  filter  of 
the  same  order.  The  Chebyshev  filter  is  an  all-pole  type  and  is  based 
on  the  nth-order  Chebyshev  polynomials,  Cn(u>) • 


cos  (n  cos”1  uj)  ,  0  <  u>  <  1  , 
cosh  (n  cosh”1  j)  ,  u>  >  1  . 


The  recursion  formula  for  finding  the  nth-order  polynomial  is 

CQ  (oi )  =  1  , 

C^cj)  =  <u  , 

C  («)  =  2u>C  .Uo)  -  C  „(w)  . 
n  n-1  n-2 

The  equations  for  calculating  the  pole  locations  are  given  in  appen¬ 
dix  A. 


2.4.3  Elliptic 

The  third  type  of  filter  approximation  considered  in  this 
report  is  called  the  elliptic  approximation  and  was  first  introduced  by 
Wilhelm  Cauer.  This  approximation  is  based  on  the  J6cobi  elliptic  sine 
functions.  In  this  approximation,  the  error  in  the  passband  is  again 
distributed  evenly  in  an  oscillating  manner.  However,  instead  of  a 
monotonically  decreasing  characteristic  in  the  stop  band,  the  stop-band 
attenuation  oscillates  between  infinity  and  a  prescribed  maximum.  Thus, 
there  is  an  equiripple  characteristic  in  both  passband  and  stop  band. 
Hie  elliptic  approximation  is  more  efficient  than  the  Butterworth  and 
Chebyshev  approximations  in  that  the  transition  between  passband  and 
stop  band  is  steeper  for  a  given  filter  order.  An  elliptic  filter 
transfer  function  has  both  poles  and  zeros  and  has  been  shown  to  be 
optimal  in  the  sense  of  having  the  sharpest  transition  of  any  approxi¬ 
mation.  The  technique  used  to  calculate  the  coefficients  for  the  ellip¬ 
tic  filter  is  given  in  appendix  A  and  was  taken  from  the  development  by 
Antoniou.1 


1A.  Antoniou,  Digital  Filters:  Analysis  and  Design,  McGraw-Hill  Book 
Co.,  Inc.,  New  York  (1979). 
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2.5  Digital  Filter  Realization 


There  are  two  possible  canonical  forms  of  linear  difference 
equations  that  can  be  realized  from  the  z  domain  transfer  function, 
H(z).  One,  the  cascade  form,  follows  from  the  factored  form  of  H(z). 
The  other,  the  parallel  form,  requires  the  expansion  of  H(z)  into 
partial  fractions.  The  parallel  form  was  chosen  for  this  project 
because  of  its  reduced  sensitivity  to  noise. 

Realization  of  a  recursive  digital  filter  in  parallel  form  is 
shown  in  figure  2. 


INPUT 


OUTPUT 


ZAnk  +  Ait,  z"1 

- ss — i; - 

K=  1  1  +  ®1k  r*  +  r* 


Figure  2.  Recursive  digital  filter  in  parallel  form. 


2.6  Band  Transformations 


For  this  project,  high-pass,  band-pass,  and  band-stop  filters 
are  obtained  from  the  digital  low-pass  prototype  by  a  transformation  in 
the  z  domain.  Another  approach  is  to  design  the  analog  low-pass  proto¬ 
type,  perform  the  transformation  in  the  s  domain  to  high  pass,  band 
pass,  or  band  stop,  and  then  transform  the  resulting  transfer  function 
to  the  z  domain  by  the  bilinear  transformation. 

The  work  of  Constantinides5  was  used  to  develop  the  z  domain 
transformations.  The  equations  are  given  in  appendix  A. 


5A.  G.  Constantinides,  Spectral  Transformations  for  Digital  Filters, 
Proc.  IEE,  U7_  (August  1970),  1585-1590. 


3.  FORTRAN  IV  PROGRAM 

3. 1  Outl ine 

A  FORTRAN  IV  computer  program  was  written  using  the  equations 
and  the  techniques  discussed  in  section  2.  The  program  was  designed  in 
a  modular  form  to  allow  for  easy  modification  as  required.  Hie  flow 
chart  of  the  program  is  given  in  figure  3.  The  main  part  of  the  program 
is  simply  a  controller  that  calls  in  subroutines  as  required.  The 
blocks  in  a  column  in  the  center  of  the  flow  chart  are  all  in  the  main 
program,  and  the  blocks  to  the  right  or  the  left  of  center  are  separate 
subroutines.  For  instance,  if  a  different  form  of  input  data  is  re¬ 
quired  or  a  different  output  format  is  desired,  then  only  the  subroutine 
dealing  specifically  with  that  function  need  be  modified.  If  a  differ¬ 
ent  analog  design  type  or  technique  is  to  be  used,  it  can  be  either 
substituted  for  one  of  the  existing  three  design  types  or  added  on  by 
extending  a  "computed  go  to”  statement  in  the  main  program. 

All  data  are  stored  in  labeled  and  unlabeled  common  areas  of  a 
subroutine  simply  by  including  the  proper  common  statement. 

Hiis  program  was  developed  by  using  an  IBM  System/370  com¬ 
puter.  The  H  Extended  Plus  FORTRAN  compiler  was  used  with  the  auto¬ 
double  option  specified.  This  option  doubles  the  overall  precision 
specified  in  the  program.  Maximum  possible  precision  is  used  for  all 
the  design  calculations  as  well  as  in  the  implementation  itself.  This 
precision  is  16  bytes  for  real  variables  and  32  bytes  for  complex  varia¬ 
bles.  The  complete  program  listing  is  given  in  appendix  B. 

The  flow  of  the  program  is  straightforward  because  it  begins  at 
the  top  and  progresses  without  any  diversions  to  the  bottom  of  the 
chart.  The  basic  steps  of  the  program  and  the  general  function  of  each 
subroutine  are  as  follows. 

3.1.1  RDATA 

After  the  program  is  initiated,  it  immediately  calls  the 
RDATA  subroutine.  This  subroutine  reads  in  the  necessary  filter  design 
criteria  and  the  data  samples  to  be  processed.  Only  batch  processing  of 
data  is  possible  with  the  existing  program,  but  the  conversion  to  real 
time  processing  would  not  be  difficult.  The  input  data  array  is  dimen¬ 
sioned  for  up  to  1000  amplitude  points  at  the  specified  sampling  fre¬ 
quency.  If  only  a  filter  design  is  required,  then  the  number  of  data 
points  (NTIMES)  is  set  equal  to  zero,  and  RDATA  bypasses  the  input  of 
data  to  be  processed.  Setting  NTIMES  equal  to  zero  bypasses  also  the 
implementation  subroutine.  RDATA  also  prints  out  the  design  criteria 
data  for  future  reference. 

Hie  specific  format  required  for  data  input  is  given  in 
section  3.2. 
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3.1.2  WARP 


The  RDATA  subroutine  returns  control  to  the  main  proqram, 
which  calculates  the  prewarped  analog  design  frequency  (WARP)  and  then 
calls  the  appropriate  analog  design  subroutine.  The  analog  filter 
subroutines  calculate  the  poles  and  the  zeros  (  if  any)  of  the  squared 
magnitude  transfer  function  that  ' ie  in  the  left  half  of  the  complex  s 
plane.  The  resulting  transfer  function  a  t  :r  t  the  right-half  plane  poles 
and  zeros  are  eliminated  is  of  the  f*  -rat 

( AMP)  (s  -  z^s  .  .  .  U  - 

H(  s)  - - - , 

(s  -  ?1)(s  -  P2)  •  .  (s  -  P,) 

where  z^  and  are  complex  conjugates  of  zx  and  p^  (for  all  i).  The 
equations  needed  to  calculate  these  poles  and  zeros  are  given  in  appen¬ 
dix  A  and  are  derived  as  discussed  in  section  2. 

3.1.3  FACTOR 


The  FACTOR  subroutine  is  then  called.  FACTOR  expands  the 
transfer  function  into  partial  fractions.  If  the  numbers  of  poles  and 
zeros  are  equal,  then  this  expansion  has  a  constant  gain  equal  to  the 
amplitude  multiplier  (AMP)  of  the  transfer  function.  FACTOR  then  com¬ 
bines  the  complex  conjugate  pairs  and  applies  the  bilinear  transfor¬ 
mation  , 


to  obtain  the  coefficients  for  the  parallel  second-order  sections  of  the 
digital  filter.  The  general  form  of  the  resulting  z  domain  transfer 
function  is^ 


H(z~l) 


(---)  f ,  *oi  *  _2 . 

(—>  1  +  BUZ_1  +  B2iZ  ^ 

i=1 


where  r  =  n/2  for  n  even,  r  =  (n  +  1)/2  for  n  odd,  and  n  is  the  order  of 
the  filter.  The  coefficients  A^,  A^,  B1if  and  B2i  are  calculated  by 
FACTOR  using  the  analog  poles  ana  zeros.  If  the  order  of  the  filter  is 
odd,  then  there  is  one  real  pole,  so  the  A^  and  B2  coefficients  are 
equal  to  zero  for  that  pole. 

^G.  C.  Temes  and  S.  K.  Mitra ,  Modern  Filter  Theory  and  Design,  John 
Wiley  &  Sons,  Inc.,  New  York  (1973). 
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3.1.4 


F  PLOT 


A  plot  of  the  frequency  res|)onsp  (magnitude  and  phase)  of 
the  low-pass  diqital  filter  is  then  made  by  the  FPLOT  subroutine.  It  is 
plotted  by  letting 

_i  -jiitT  m  m 

z  1  =  e  -  cos  (i)T  -  •)  sin  u>T 

in  the  z  domain  transfer  function  given  above. 

The  plots  given  in  this  report  were  made  on  a  Houston 
Instruments  plotter  using  existing  plotting  software.'’  Both  the 
amplitude  and  phase  responses  are  plotted  over  the  normalized  range  of 
frequencies  of  zero  to  one.  That  is,  the  frequency  scale  is  normalized 
to  the  sampling  frequency.  FPLOT  also  prints  the  values  of  the  digital 
filter  coefficients  to  25  decimal  places  for  ’uture  use. 

3.1.5  Filter  Choice 


The  main  proqram  then  determines  if  a  high-pass,  band-pass, 
or  band-stop  filter  is  desired  or  if  the  low-pass  filter  is  ready  to  be 
implemented . 

The  transfer  function  for  the  required  high-pass,  band-pass, 
or  band-stop  filter  is  obtained  by  r-  .ppropriate  transformation  in  the 
z  domain  of  the  low-pass  digital  filter.  Band-pass  and  band-stop  trans¬ 
formations  result  in  a  doublinq  of  the  filter  order. 

After  tr  u.sformation ,  the  frequency  characteristics  of  the 
new  filter  are  plotted  by  FPLOT. 

3.1.6  FILIMP 


The  actual  implementation  of  the  digital  filter  is  simple. 
The  FILIMP  subroutine  performs  this  function  and  processes  the  data  to 
be  filtered.  Application  of  the  inverse  z  transform  to  the  general  form 
for  the  z  domain  transfer  function  results  in  the  recursive  difference 
equation 

Yout(kT)  =  X(kT)  +  X{(k  -  1)T1  , 


' Thomas  V.  Noon  and  Egon  Marx,  User's  Manual  for  the  Modular  Analysis- 
Package  Libraries  ANAPAC  and  TRANL,  Harry  Diamond  Laboratories  HDL-TR- 
1782 -S  (September  1978). 
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where 


r 

X(kT)  =  l  Xi(kT)  , 
i=1 

Xi(kT)  =  AoiYin(kT)  +  AuYin[(k  -  1  )T] 

-  BuXi[(k  -  1)T]  -  B2iXi[(k  -  2  )T]  . 

3.1.7  PDATA 

The  PDATA  subroutine  provides  plots  of  the  input  and  output 
data.  This  routine  may  be  easily  modified  to  provide  the  output  data  in 
any  desired  format. 

3.2  Program  Use 

There  are  at  most  nine  parameters  that  must  be  supplied  to  the 
program  for  the  design  of  a  filter.  The  minimum  number  needed  is  five 
for  the  design  of  a  Butterworth  low-pass  filter.  This  minimum  set, 
which  is  used  in  all  designs,  consists  of  these  parameters: 

a.  ITYPE ,  the  type  of  analog  design 

1  =  Butterworth 

2  =  Chebyshev 

3  =  Elliptic  (Cauer) 

b.  N,  the  order  of  the  filter 

c.  ITRANS,  corresponding  to  the  band  transformation  desired 
0  =  None 

1  =  High-pass  transformation 

2  =  Band-pass  transformation 

3  =  Band-stop  transformation 

d.  FC,  the  cutoff  frequency  desired  ( in  hertz) 

e.  FS,  the  sampling  frequency  of  the  data  (in  hertz)  to 

which  the  filter  will  be  normalized 
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If  a  low-pass  Chebyshev  is  required,  then  an  additional  param¬ 
eter,  EPS  1 ,  the  minimum  allowed  amplitude  in  the  passband,  also  must  be 
entered.  Similarly,  for  elliptic  filters,  EPS  1  must  be  given  along  with 
EPS2,  the  transition  region  selectivity  factor.  The  maximum  value  for 
EPS2  is  0.95.  All  of  these  parameters  are  read  from  the  same  data 
card.  These  are  the  formats  for  all  the  data  cards  and  the  order  in 
which  the  cards  are  read  by  the  program : 

ITY PE , N , ITRANS , FC , FS , E  PS  1 , EPS  2 
Format — 315,  4E10.3 

If  ITRANS  =  0,  skip  reading  FI  and  F2. 

F 1 ,  F2  , 

Format — 2E10. 3 

NTIMES 
Format — 15 

If  ITIMES  =  0,  skip  reading  VIN. 

VIN 

Format — 8E10.3 — maximum  of  125  cards 
Heading  for  low-pass  design 
If  ITRANS  =0,  skip  next  heading. 

Heading  for  band  transformation 
If  NTIMES  =  0,  skip  next  two  headings. 

Heading  for  VIN  plot 
Heading  for  VOUT  plot 

If  a  band  transformation  is  specified,  then  an  additional  data 
card  is  needed  that  has  the  lower  cutoff  frequency,  FI  (also  the  high- 
pass  cutoff  frequency) ,  and  the  upper  cutoff  frequency,  F2,  for  the 
band-pass  and  band-stop  filters. 

Next,  the  number  of  data  sample  points  is  read.  If  NTIMES  is  set 
equal  to  zero,  then  no  input  data  are  read,  and  the  program  is  termi¬ 
nated  after  completion  of  the  filter  design. 

Separate  titles  are  read  in  for  the  low-pass  filter  transfer 
function  plots,  any  band  transformation  frequency  plots,  and  the  heading 
to  appear  on  the  data  plots.  These  titles  can  be  up  to  80  characters  of 
any  form  desired. 
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A  sample  filter  design  is  given  now  to  illustrate  the  order  and 
the  manner  in  which  the  data  are  assembled.  A  fourth-order  elliptic 
band-pass  filter  is  chosen  since  this  design  requires  input  of  all  the 
design  parameters.  We  have  this: 

ITYPE  =  3 

N  =4 

ITRANS  =  2 

Now  let  us  assume  that  the  sampling  frequency  is  10,000  Hz  and  that  the 
passband  is  to  be  from  1000  to  3000  Hz.  In  addition,  we  allow  the 
passband  amplitude  to  vary  between  1.0  and  0.9  and  minimize  the  transi¬ 
tion  between  the  passband  and  the  stop  band.  This  action  means  that  the 
cutoff  frequency  for  the  low-pass  design  must  be  2000  Hz,  which  is  the 
width  of  the  passband.  Therefore,  the  remaining  parameters  are  these: 

FC  =  2000. 0 

FS  =  10000.0 

EPS  1  =0.9 

EPS2  =  0.95 

FI  =  1000.0 

F2  =  3000.0 

These  data  are  given  in  figure  4  as  they  appear  on  the  input  data 
cards.  The  plot  title  cards  also  are  shown  in  the  figure.  If  only  a 
low-pass  design  is  requested,  then  the  second  title  card  is  omitted. 
Similarly,  if  no  data  points  are  given,  then  the  third  title  card  as 
well  as  the  input  data  cards  are  omitted. 

The  resulting  printout  for  this  example  problem  is  given  in 
figure  5.  The  low-pass  transfer  function  is  given  in  figures  6  and  7, 
and  the  band-pass  function  is  given  in  figures  8  and  9.  Running  this 
example  required  6.85  s  in  the  central  processing  unit  (CPU)  to  compile 
the  FORTRAN  coding,  2.44  s  in  the  CPU  to  link  and  edit  the  object  mod¬ 
ules  with  library  functions,  and  about  2  s  in  the  CPU  to  perform  the 
actual  calculations.  The  filtering  of  101  data  points  required  less 
than  2  s  of  CPU  time.  Additional  examples  of  filter  designs  are  given 
in  section  4. 
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Figure  4.  Input  data  for  example  1 . 


1.0417444015054475*315774110-0]  1  ■  *10b7l53*4j9  4l*1540005*0*0-0l  -4 .4624033  I  454 7 1 1 0 74 34440  I 05Q-01 

-4.1  994047544694  34  742  73944  14Q-02  >  .4700 70  12  1444  14  49 74  748  14 24(1-0*  -4.  24  3  704444 74  015  » 1444 1  729400-01 

IHF  •  1. 442991 99944? 795403*4! >7700-01  A?  •  7. 


1.44 1413 307774*041 74  514  36210-01 
9. 2440  «0  1U  1411  1424  7 15405 1  30-01 


-1.79 1435994 141  014 124414 7414  0*01  - 7 .49 15 1 9 1444 10 744440 11944 1 10 -02  -1 .251141477457449 14529445440*00 
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4.94444029475540514114711100.01 
4.72407451947145744554541400-01 
5.7214444  74  7120  11  7770404  4  720-01 
9. 554434 10 720 1742504 5 7224400-01 


Figure  5.  Output  data  for  example  1. 


Figure  6.  Low-pass  Cauer  filter:  amplitude  versus  frequency 
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Figure  7.  Low-pass  Cauer  filter:  phase  versus  frequency 
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4.  DESIGN  EXAMPLES 


Several  examples  are  given  in  this  section  to  illustrate  the  differ¬ 
ences  in  the  various  analog  designs  and  also  to  show  the  results  of  band 
transformations.  In  addition,  a  simple  waveform  consisting  of  three 
sinusoids  at  different  frequencies  is  used  to  demonstrate  the  filtering 
action  of  the  various  digital  filters.  . 

4. 1  Analog  Design 

The  three  types  of  analog  designs  are  discussed  in  section  2. 
A  graphical  presentation  of  the  basic  differences  in  these  three  designs 
is  given  here.  A  fourth-order  elliptic  low-pass  design  is  given  in 
section  3.2.  Similar  fourth-order  digital  filters  using  the  Butterworth 
and  Chebyshev  analog  designs  are  presented  in  figures  10  to  13.  Compar¬ 
ing  these  transfer  functions,  we  can  easily  see  that  the  Butterworth 
filter  offers  the  smoothest  response  in  both  the  passband  and  the  stop 
band.  The  sharper  transition  from  passband  to  stop  band  of  the 
Chebyshev  design  over  the  Butterworth  is  obtained  at  the  expense  of 
smoothness  in  the  passband.  The  elliptic  design  obviously  has  the 
sharpest  transition  of  the  three  designs.  However,  the  ripple  in  both 
the  passband  and  the  stop  band  is  the  penalty  for  this  sharp  transi¬ 
tion.  The  design  to  use  is  chosen  according  to  the  requirements  of  the 
specific  filtering  problem  being  considered. 

4. 2  Band  Transformations 

An  example  of  a  band-pass  transformation  is  given  in  figures  8 
and  9.  The  low-pass  Butterworth  filter  of  figures  10  and  11  was  trans¬ 
formed  to  a  high-pass  filter  with  a  cutoff  of  3000  Hz  as  shown  in  fig¬ 
ures  14  and  15.  The  band-stop  dual  of  the  passband  filter  was  made  by 
using  the  Chebyshev  design  and  is  shown  in  figures  16  and  17. 
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Figure  14.  High-pass  Butterworth  filter:  amplitude  versus  frequency 
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Figure  15.  High-pass  Butterworth  filter:  phase  versus  frequency 
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Figure  16.  Band-stop  Chebyshev  filter:  amplitude  versus  frequency 
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Figure  17.  Band-stop  Chebyshev  filter: 

phase  versus  frequency 

4. 3  Data  Processing 


A  simple  waveform  was  generated  to  illustrate  the  operation  of 
the  filters  on  input  data.  The  waveform  used  was  the  sum  of  three 

sinusoids  at  0.05,  0.15,  and  0.35  times  the  sampling  frequency  with  peak 
amplitudes  of  3,  2,  and  5,  respectively.  This  input  waveform  is  shown 

in  figure  18.  A  20th-order  Butterworth  low-pass  design  was  used  to 
eliminate  the  two  highest  frequency  sine  waves,  as  shown  in  figure  19. 
A  filter  of  high  order  was  used  to  demonstrate  the  group  delay  effect  of 
the  filter.  This  effect  shows  at  the  beginning  of  the  output  waveform 
as  a  delay  time  before  any  output  occurs.  Next,  a  fourth-order  Butter- 
worth  high-pass  filter  was  used  to  eliminate  the  two  low-frequency 

signals.  This  result  is  shown  in  figure  20.  The  low  sampling  rate 

caused  the  distortion  on  the  sine  wave  seen  in  this  figure.  The  initial 
design  example  was  used  to  pass  only  the  center  sinusoid  as  shown  in 
figure  21.  Here,  again,  the  effect  of  the  low  sampling  rate  is  seen  as 
a  distortion  of  the  waveform.  However,  the  effect  is  not  as  great 

because  of  the  lower  frequency  of  the  sinusoid.  The  Chebyshev  band- stop 
dual  of  this  passband  filter  was  then  used  to  eliminate  the  middle  sine 
wave.  This  waveform  is  given  in  figure  22. 
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Figure  18.  Input  waveform  obtained  by  summing  three  sinusoids. 
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Figure  19.  Output  waveform  obtained  by  using  low-pass  Butterworth 
filter . 


Figure  20 . 


Output  waveform  obtained  by  using  high-pass  Butterworth 
filter. 
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Figure  21.  Output  waveform  obtained  by  using  band-pass  Cauer  filter 
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Figure  22.  Output  waveform  obtained  by  using  band-stop  Chebyshev 
filter. 


5.  CONCLUSIONS  AND  RECOMMENDATIONS 


A  versatile  program  has  been  developed  for  designing  and  implemen¬ 
ting  recursive  digital  filters.  Three  well-known  analog  filter  designs 
were  used  with  the  bilinear  transformation  to  obtain  the  desired  digital 
filters.  The  bilinear  technique  was  chosen  because  it  eliminates  the 
effects  of  aliasing  that  occur  with  other  techniques.  However,  the 
resulting  warping  of  the  frequency  scale  may  not  be  acceptable  in  some 
applications.  The  extent  to  which  this  warping  affects  the  resulting 
output  was  not  investigated.  Examples  using  the  different  analog  de¬ 
signs  are  given  in  section  4  along  with  transformations  from  low-pass  to 
high-pass,  band-pass,  and  band-stop  designs.  All  of  the  filter  design 
calculations  as  well  as  the  implementations  use  the  highest  precision 
possible  for  the  IBM  System/370  computer.  No  investigation  was  made  to 
determine  the  effects  of  using  lower  precision. 

All  of  the  equations  used  in  developing  the  code  are  included  in 
appendix  A,  and  the  complete  computer  program  code  is  listed  in  appendix 
B.  A  bibliography  of  related  publications  also  is  included.  This 
bibliography  is  by  no  means  complete  since  the  amount  of  information 
published  seems  almost  endless.  Although  there  is  a  large  amount  of 
information  available,  no  one  source  proved  adequate  in  presenting  all 
of  the  steps  necessary  to  completely  develop  a  digital  filter.  It  is 
hoped,  therefore,  that  this  document  will  prove  helpful  to  others  ven¬ 
turing  into  this  field  for  the  first  time. 

Although  the  program  developed  here  is  completely  operable,  there 
are  many  ways  in  which  the  program  may  be  improved  or  modified  to  better 
suit  a  particular  filtering  need.  Also,  investigation  and  analysis  of 
the  various  errors  associated  with  the  techniques  used  should  be  con¬ 
ducted  . 
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APPENDIX  A 


The  following  equations  were  used  in  the  development  of  the  com¬ 
puter  code  for  designing  and  implementing  digital  filters.  All  of  the 
equations  are  referenced  in  the  main  body  of  this  report.  The  applica¬ 
tion  of  each  set  of  equations  is  given  in  the  flow  chart  of  the  program 
(fig.  3). 

A-1 .  PREWARPING  (WARP) 


A- 2.  8UTTERW0RTH 

The  transfer  function  for  a  Butterworth  low-pass  filter  of  order  n 
is  given  by 


where 


H_ 


H(s)  = 


0 


",  (s  - 


H„  =  ui  , 
0  c 


=  cutoff  frequency, 

~  j[(*/2)+u(2k-1)/2nl  ,  ,  . 

Pk  =  uce  '  k  =  2' 


n . 


A- 3.  CHEBYSHEV 

The  squared  magnitude  function  is  of  the  form 


|H( ju) |2 


1 

1  +  E2C2(w) 


where 


e  =  ripple  factor, 

Cn ( (u )  =  Chebyshev  polynomial. 
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ELLIPTIC  (CAUER) 


For  the  normalized  filter  u)c  =  1, 
n  =  filter  order. 


k  =  selectivity  factor  =  “p/wa, 

Ap  =  maximum  passband  attenuation  (in  decibels), 
Ag  =  minimum  stop-band  attenuation  (in  decibels). 


Of  the  parameters  Ap,  Aa,  k,  and  n,  only  three  are  independent. 
If  three  of  the  four  are  specified,  the  fourth  is  automatically  fixed. 
For  this  project,  n,  k,  and  Ap  are  specified: 

k  =  EPS  2, 

Ap  =  20  log  ( 1  -  EPS1) . 


Aa  is  then  given  by 


A  =10  log 

3 


-0.  1 A 

10  P  -  1 


+  1 


16qr 
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H(s) 


\ 


\ ( n+1 ) /2  ,  n  odd  , 
I n/2 ,  n  even  , 


s  =  2  I 
S  T 


-  z  ' 


1  +  z 


H(z~l)  =  (l  +  z~l)  ^ 


'Ok 


+  A1kz 


-1 


1  +■  B,.z'*  +  B„,  z 


k=1 


Ik 


2k 


pk = 

_Jk  + 

k  ' 

A0k  = 

M'  * 

V2)  -  MV2)]20 

A1k  = 

-  DvO 

-  V2)  *  UV2)] 

B1k 

-2[’  -( 

V2)2  -  (V2)2]'0 

B2k  = 

[(’  -  «k 

/2)2  *  (V2)2]/D  - 

D  -  0 +  v2)2  +  (v2)2  • 


A-6.  MAGNITUDE  AND  PHASE  ( F PLOT ) 


The  frequency  response  of  the  designed  digital  filter  is  calcula¬ 
ted  by  substituting  e  -'uiT  for  z  *  in  the  transfer  function  and  then 
using  complex  arithmetic  to  obtain  the  magnitude  and  the  phase  as  a 
function  of  frequency. 


s"1  -  e"3wT  = 


W  =  cos  w T  -  j  sin  wT  , 
r 


k=1 


V  A0k  +  A1kW 
H(W)  =  (1  +  W)  )  -  . 

^  1  *  S1kW  +  B2k«2 
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HIGH-PASS  TRANSFORMATION 


APPENDIX  A 


A  low-pass  digital  filter  of  cutoff  frequency  0  is  first  designed 
for  the  desired  order,  n,  and  type.  Let  ojp  =  the  desired  high-pass 
cutoff  frequency. 

Then  substituting  for  z  ^  in  the  low-pass  transfer  function, 

,  z  '  +  a 


1  +  uz‘ 


where 


(K±) 

(^) 


results  in  the  transfer  function  of  the  desired  high-pass  filter. 


Z  Ok 


A0k  +  A1kz  1 


+  B2kz  " 


A0k  =  (1  '  «> 


A0k  aA  Ik 


A1k  =  51  "  a) 


afl0k  "  A1k 


B1k  =  [2«  "  B1k(1  +  a2' 


+  2B2kaJ/D'  , 


=  (a2  - 


aB1k  +  B2k 


0'  =  1  -  aB1k  +  a2B2k  . 


A-8.  BAND-PASS  TRANSFORMATION 


A  low-pass  digital  filter  of  cutoff  frequency  6  is  first  designed 
for  the  desired  order  and  type.  The  cutoff  frequency  equals  the  band¬ 
width  of  the  desired  band-pass  filter,  F2  -  FI,  where  F2  and  FI  are  the 
upper  and  lower  cutoff  frequencies  (in  hertz),  respectively. 
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A  substitution  for  z-^  in  the  low-pass  transfer  function  results 
in  the  transfer  function  of  the  desired  band-pass  filter. 

In  general,  this  transformation  is1'2 


>-l 


2uk 
k  +  1 


k  -  1 
k  +  1 


k  +  1 


1  z"2 


2ak 
k  +  1 


+  1 


where 

COSH^F1  +  F2)T 
u  _  COSIT  ^F2  -  F^T' 

k  =  cot  [h(F2  -  Fi)t]  tan  (nfJT). 

By  letting  k  =  1,  then  &  =  F2  -  FI,  and  the  transformation  is  simpli¬ 
fied. 


z 


-1 


The  resulting  band-pass  transfer  function  has  the  form 


n/2 


A11kz  1  +  A01k 


C1  /  UK 


B2 1kz 


-2 


A12kz  1  “  A02k 
1  —  B ^ 2)^ z  ^  ^  ^ 


22k 


*A.  G.  Constant inides ,  Spectral  Transformations  for  Digital  Filters, 
Proc.  IEE,  117  (August  1970),  1585-1590. 

2A.  Antoniou,  Digital  Filters:  Analysis  and  Design,  McGraw-Hill  Book 
Co.,  Inc.,  New  York  (1979). 
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for  n  even.  For  n  odd,  there  is  one  real  pole,  and  the  transfer  func¬ 
tion  is  of  the  form 


A0[1  -  z 


-2  ' 


1  -  a[  1  -  BtJ z"1  -  B i : 


-2 


(n-1 )/2 


r-'  /  A11kz  1  +  A01k 

k=1  V  +  1  +  ®21k2 


(1 


A1 2kz  1  "  A02k 
1  “  B12kz  1  +  B22kz 


A-9.  BAND-STOP  TRANSFORMATION 

The  procedure  for  the  band-stop  transformation  is  the  same  as  for 
the  band-pass  transformation.  In  general,  the  transformation  is 


where 


,-l 


2a  1  -  k 

-2  -  Z  1  +  - 

1  +  k  1  +  k 


1  -  k  _  2a 

-  z  ^  -  - 

1  +  k  1  +  k 


+  1 


COSH^F1  +  F  \T 

U  =  - - - , 

cosw^F2  -  F^T 

k  =  tan  [h(f2  -  F^t]  tan  (tt3T). 


By  letting  k  =  t,  then  6  =  Fg/2  -  (Fj  -  F.j),  and  the  transforma¬ 
tion  is  simplified. 
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The  resulting  band-stop  transfer  function  has  the  form 


hbs(z_1)  =  l1  -  2az 


-1 


+  z 


-2 


^  /  A11kz  1  "  A01k 

^  V  B11kz  1  + 


2  Ik 


,~2 


A12kz  1  +  A02k 
1  +  B12kz  1  +  B22kz  ^ 


for  n  even.  For  n  odd,  the  transfer  function  has  the  form 


hbs(z_1) 


( 1  -  2az-l  +  z“2 
1  -  a(  1  +  Bjz-1  +  B.,z-2 


+  (  1  -  2az  1  +  z” 


(n-1 }/2  /  -I  .  j' 

Zf  A11kz  A0 Ik 

K=)  V  -  B11kz_1  +  B2  Ik2 


A1 2kz  1  +  A02k 
1  +  B1 2k z  1  +  B22kz  2' 


A- 10.  FILTER  IMPLEMENTATION 

From  the  z  domain  general  transfer  equation 


V-1  A0i  +  A1iz  1 

H(z'l)  =  (1  +  z-1)  )  - - - - 

Z_ I  1  +  B..Z  ^  +  B,.z  L- 

i=1  11 


a  time-domain  difference  equation  is  obtained.  The  difference  equation 
is  then  used  to  process  data  samples. 

Yout(kT)  =  X(kT)  +  xt(k  -  1 )T]  , 
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where 


X(kT) 


Xi(kT) 


l  X^kT), 
i=  1 


AY.  (kT)  +  A.. Y,  [(k  -  1 )T] 
Ol  xn  lx  in 


-BuXi[(k  -  1  )T]  -  B2±Xi  [  ( k  -  2)T]  , 
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APPENDIX  B. — FORTRAN  IV  COMPUTER  PROGRAM  FOR  RECURSIVE  DIGITAL  FILTERS 
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FriSCiDLNG  PiQ*  BUMK-NOT  Fli>SD 


CONTENTS 


MAIN  PROGRAM  .... 
SUBROUTINE  RDATA 
SUBROUTINE  BUTTER 
SUBROUTINE  FACTOR 
SUBROUTINE  CHEB  . 
SUBROUTINE  CAUER 
SUBROUTINE  HIGHP 
SUBROUTINE  FPLOT 
SUBROUTINE  BAN DP 
SUBROUTINE  BANDS 
SUBROUTINE  FILIMP 
SUBROUTINE  PLT DAT 
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The  main  program  and  all  of  the  subroutines  for  recursive  digital 
filters  are  included  in  this  appendix.  The  only  exceptions  are  standard 
library  functions  and  the  plotting  routine.  If  a  machine  other  than  an 
IBM  is  used,  it  may  be  necessary  to  modify  the  format  and  data  state¬ 
ments  . 


B— 1 •  MAIN  PROGRAM 


REALMS  40*20),  41*20),  B  1*20),  B2*20)t  PI ,A2 
REAl«8  AMP 

CQMPLE  X* 16  P  *2D ) ,  l *20) 

COMMON  FC,  FS,  EPSl,  FPS2,  PI,  1TYPE,  N,  ) TRANS 
C0MM0N/VDATA/V1N*1000I,  V0UT  flOOO)  ,NT IMES 
CQHM0N/TRANS/  FI ,F2 

COMM ON /F ILT /WC,  N1,W2,AMP,A0,A1,A2,  Bl,  B2,  P,  I 
C  ITVPE-  FILTER-  BUTTER  WORTH  *  1 ) ,  CHEB YSHE VI2 ) ,  ELLIPT1C*3) 

C  N  -  ORDER  OF  FILTER,  20  MAX. 

C  1  TRANS-  NON  E  *  C)  ,  dPIll,  BP*2),  BSI3) 

C  FC-  COTOFF  FREQUENCY 

C  FS-  SAMPLE  FREQUENCY 

C  EPS1  -  MINIMUM  4MPLITU0E  IN  PASSBAND 

C  EPS2  -  TRANSITION  COEFFICIENT  FOR  CAUER,  0.95  MAX. 

C  FI  -  LONER  CUTOFF  FREQUENCY,  PASS  AND  STOP  FILTERS  ONLY 

C  F 2  —  UPPER  CUTOFF  FREQUENCY,  PASS  AND  STOP  FILTERS  ONLY 

C  NT IMES  -  NUMBER  OF  DATA  POINTS 

C  VIN-  INPUT  DATA 

C  VQUT-  OUTPUT  DATA 

P)s4.DDO*DATANtl.DDO) 

CALL  RDATA 
MC=2*DTAN*PI *FC/FS» 

M 1 =P I *F 1 /FS 
N2-PHF2/FS 

GO  TO  *10,20.301, ITYPE 
ID  CALL  BUTTER 

tO  TO  40 
2D  CALL  CHEB 

CO  TO  40 

3D  CALL  CAUER 

4D  CALL  FACTOR 

CALL  FPLOT 
4J-ITRANSM 

GO  TO  1100,50,60,70 !,JJ 
5D  CALL  HIGHP 

GO  TO  60 
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B— 1 .  MAIN  PROGRAM  (Cont'd) 


to  CALL  BANDP 

GO  TO  83 

73  CALL  BANOS 

83  CALL  F  PLOT 

133  lFiNTlHES.EO.OlG3  TO  110 
CALL  F  11  1HP 
CALL  P  L  T  DAT 
110  CONTINUE 

STOP 
END 


B-2.  SUBROUTINE  RDATA 


SUBROUTINE  RDATA 
REAL#8  PI 

COMMON  FC,  FS,  EPSI  •  EPS?*  PI,  ITYPE,  N,  1  TRANS 
COMNON/VDATA/vmiOOO  1,  VOUT 11000 1 ,NT I  ME S 
COMMON/TRANS/  F  1  ,  r  2 

REAOiS.lOl  HYPE,  N,  ITRANS,  FC.FS.EPSl.  EPS? 

13  FORMAT (315*  4E13.31 

F  2-  0 
F  1  =F  2 

IF  (ITRANS  .  E  3  .0 1  GO  TO  100 
READ«5.?0»  F1.F2 
23  FORMAT  t?E10 .3) 

READI5,25)NTIMES 
2S  FORMAT  1 1  SI 

IF  INTJMES. Efa.OIGD  IQ  50 
JJ=  tNT IMES«7)/8 
DO  100  K -1 »  J J 
L« IK-1 )»8«1 
M*K  *8 

130  READfS*301IVINII 1*1 =L»MI 

33  F  ORMAT I8E10.31 

53  CONTINUE 

NR1 TEIb,45  1 

NS  FORMAT  l////,5X5tmYPE  *7X  1HN  ,6XfcH  I  TRANS.9X2HFC  *  13X2HFS  ,  1 2X4HEP  S 1 , 
111X4  HE  PS2*12X?MF1*1  3X2HF21 

UR1TE(6,40)  I  TYPE •  N,  ITRANS*  FC ,  FS,  EPSI.  EPS2 »  FI,  12 
43  FORMATI/, 3I5XIS1*  b I5XE 10.31  ,//// 1 

RETURN 
END 
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B-3.  SUBROUTINE  BUTTER 


SUBROUTINE  BUTTER 

R  E  A  L  *8  A  0  4  2  0  1  ,  41420),  BII20),  B2420),  Pl.PrtI 
C3NPLEX*16  P  4  20J,  2  4201 

COMMON  EC,  F  S  »  EPS1,  EPS2,  PI,  IT  YPE  ,  N,  I  TRANS 
COMMON/F 1LT/WC,  U 1 , U2  ,  AMP, AO  ,A 1 ,  A2  ,  SI,  82,  P,  l 
REAL*8  A  MP , A  2 
A*N 

K  =  N~2X INTIN/2.1 
DO  20  1*1, N 
B  =  I 

1FIK.EQ.0I  CO  TO  S 

PHI -PI *4 (A«l )/2.j«B-l -Ol/A 

GO  TO  ID 

PHl=Pl*41.U/12.0*A)*fB-l.Q)/A«0.5I 
Pill-  WC*OCMPIXIOCOSI PHI I.OSINIPHI )) 

CONTINUE 
AMP  =WC  **N 
RETURN 
END 


B-4.  SUBROUTINE  FACTOR 


SUBROUTINE  FACTOR 
RE  AC *8  AMP 

REAl*8  A0I2  ),  411201,  BII20),  B2I201,  PI,  RR , R 1 ,0 ,  AK,  BK,A2 

COMPLE  X*lb  P I  20 1  ,  2  420),  R420),  X 

COMMON  FC,  FS,  EPS1,  EPS2,  PI,  1  TYPE ,  N,  ITRANS 

COMMON/F ILT/KC,  N 1 , M2  ,AMP, AO ,A 1 , A2  ,  81,  B2,  P,  2 

A2=1.0D0 

K*  I  NT  4IN*l)/2.) 

DO  50  1*1, K 

X=OCMPLX(1„000,  0.0001 
00  20  Jx 1 , N 
IF  4J.EQ.1)  SO  TO  20 
X  =  X*IP1I)-PI  Jl) 

20  CONTINUE 

R ( 1 1  -  AMP/X 
IF4N.E0.11G0  TO  30 
X  =  OCMPLX  41 .ODO.D.ODOl 
1 F  4  I  TYPE .NE . 3  1G3  TO  50 
KK  =  2*fN/2» 

00  30  J=1,KK 
X=X*|P|I  l-Zt  Jl) 
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B-4.  SUBROUTINE  FACTOR  (Cont'd) 


3D  CONTINUE 

Rf !)=RII  )*X 
AMP^AMPMKK  +  l-N) 

53  CONTINUE 

I F  C  I  TYPE  .NE  .  3  )A1P=0  .000 
00  100  1  =1,  K 
RR=DREAL  (MI)) 

RI =0  JMAG  IRI  I  I  ) 

AK  =  -DREALJP  I  I  ))*0.5DO 
BK  =  01  MAGI P  t I )) *3 . 5  DO 
0=11. 003*  AK  )**2*BK **2 
AOI I )= IRR«I 1  .0D3*AK  )  -RI*BK!/D 
Al( I ) =- I RR* ll.33D-AK)*Rl *BK )/0 
BM  1  )=  -2.000*11  .  DDQ— AK**2— BK**2 I/O 
B2U)=  ((  1.0DD-AA)**2«BK**2»/D 
130  CONTINUE 
I =2*K-N 

1FI1.E0.0IG0  TO  113 
A0IK)=A3(K)/2 
Blf K)=B1 (K)/2 
A 1  (  K  )  =  0  • 

B2 ( K ) =  3 . 

110  RETURN 
END 


B-5.  SUBROUTINE  CHEB 


SUBROUTINE  CHEB 

COMM ON /F I L  T /DC,  Ml , W2  , AMP, AO ,A1  ,A2  ,  Bl,  B 2,  P,  l 
COMMON  ft,  FS,  EPS),  EPS2,  Pi,  I  TYPE ,  N ,  1  TRANS 
C0MPIEX*16  M20),  ZI20) 

RE  AL  *8  A  0 1 2  L  )  •  AM20),  B1I20I,  B2(201 
REAL *8  CHV,SHV,X,J,AMP,Pf,A2 
EPS1=S0RT(1  .  3/EPS1* •2-1.0) 

A2*l  .003 
A  =  N 

V  =  OLOGtl .0 00/ EPS1*0 SORT  11.000*1 .ODO/E  PS1**2 ) ) /A 

CHV  =  (DEXPIV  )  *DEXP»-VM/2  .0  00 

SHV MDEXPtV  )-OEXP  I -V)  1/2.000 

JJ=»N*l)/2 

DO  50  1=1, N 

U=P1»(2.000*I1*N)-1 .000)/! 2.0D0*A) 

PI)  )  =riC*0CMPLX(0SIN (U)*SHV,  -DC DSIU ) *CHV 1 
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B-5 .  SUBROUTINE  CHEB  (Cont'd) 


53  CONTINUE 

K  *  N  —2*  JJ 
AMP  =  1.000 
If (K.EQ.OIGL  TO  93 
JJ=JJ-1 

IF (N.EQ. 1 1GO  TO  S3 
DO  75  1*1, JJ 

AHP=AMP*C0REAICP(II  I • *2* D| MAG  I P 1 1 > >**2 ) 
75  CONTINUE 

83  AMP*AMP»0REALIPI JJ*1I> 

GO  TO  125 

93  DO  100  1*1,  JJ 

AMP  =  AMP -MORE  All  Pill  I  **2«  DIMAGl  P 1 1 )  >*®2> 
130  CONTINUE 

AMP  =  AMP/DSGR T (1 .0D0*EP51  *»2> 

125  RETURN 

END 


B-6.  SUBROUTINE  CAUER 


SUBROUTINE  CAUER 

COMMON/f 1LT/MC,  Ml , W2 , AMP, AO ,A 1 ,A2 ,  Bl,  B 2,  P,  7 
COMMON  EC,  fS,  EPS1,  EPS2,  PI,  1  TYPE ,  M ,  1  TRANS 
REALMS  AMP,C*A,S0»N  ,0  ,V 

REALMS  A0(2G1,  A1I20),  811201,  B2I20),  P1,A2 
COMPLE  X*  16  2  120), PI20) 

SQ  =  0  .000 
M  *0.000 

G*C1  .30a-08LE<EP52M*2M*10.25» 

C  =  0.5DO*m.ODO-Q)/  n  ,0D0»8)) 

AMP  =  1 . 000 

G* Q«2.0*Q**5*  15.0  *Q  **9-«l  50  .0 3 
A  =  08LEIEPS1  ) 

A *010G ( f 1 .000/A* 1 .0  DO )/( 1.000/ A-l .ODO ) ) / (2 .0 *N ) 
0*0.000 
OQ  10  1*1,100 
M  *  1  — 1 

S0*S0*t-1.0D0l**M*0*MM*n*DSlNHm*M)*AI 

If CN.FQ.SOJGD  TJ  ^0 

10  w*so 

23  DO  30  1*1,100 

0*0* 1-1. 000 1 •*I*g** |) *1>*DC0SH|2.0*1*AI 
I f (Q.EO.M)GO  TO  *3 
33  M*0 
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B-6.  SUBROUTINE  CAUER  (Cont'd) 


41  S0-DABSI2.  J*Q**I0 .2 5)* SO/ 11 .000*2.0*0 )) 

C  =  N 

K=DSQRT I  II  .oDO*EPS2*50** 21*11. 0D0*S0**2/EP 521) 
IFIN.E0.11G0  TQ  100 
J J  =  N/2 

U= I2*JJ* 1-N 1/2.3 
DO  100  1*1,  JJ 
0  =  0. 0D0 

V  =  0 

DO  50  J*  1  , 1  G  0  ,2 

M  =  J— 1 

0=1  -1.0)  **M*0*MM*  J  )*DS1NI  ( N* J 1 *P l * IF LO AT  II  1-01/0*0 
N  =  J*1 

U=|-1.0) **J*Q**(M*  J)*DSINI  IH*J)*P1*IFL0AT  II )— U)/C)*Q 
IF  (O.EQ.  VIGO  TO  63 

51  V  =  0 

63  A=0 .ODO 

V  =  A 

DO  70  J=  1 , 1  Q 0  ,2 

M  =  J 

V  =  V*l-1.0I**N*Q»*IH*f1  l*DC0St2.0*«*Pl*IFLDAT»I)-U)/C) 
K  =  J*I 

V=V*I-1.01**M*Q**IM*M 1«DC0S»2.0*M*P1*IFIQATI ll-Ul/Cl 
1FIV.EQ.AIGG  TO  83 
73  A  =  V 

80  0  =  2. 0*0* *10 .251* 0/1 1.000*2.0*0 

Y=ll.0O3-EP  52*0*01*  II .00 0-0*0/ EP 52 ) 

V=DSQRTIV) 

Zll  )=DCNPlXI0.030,l  .0  00/0)*ilC 
ZIN«1-1)-DC0NJGIZ(1  >1 
A=  1  .ODO*  I  SO  *  0  )**  2 
PII)=DC!iPLX|-53»V/A  ,D  *W/ A 1  *Mt 
PlN*l-l)*DCONJG{Pti )) 

AHP=AMP*0*0* l(S3*V) *«2*I0*KI «*2)/A**2 
130  CONTINUE 

K=N-2*«N/2) 

I  FIK  .EQ.OIGO  TQ  150 

PIIN*11/ 21  =  DCHPL  X I -SO  *MC  *0 .ODD  1 

AHP  =AMP*SO*NC 

GO  TO  200 

150  AMP=AMP*EPS1 
230  CONTINUE 
RETURN 
END 
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B-7.  SUBROUTINE  HIGHP 
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SUBROJT I NE  HIGH? 

RE Al*B  A 0*2  )»  H 12  0)  •  81(201,  B2I20I,  PI 
REAL*8  A  MP , A  2 

RE  A  L  *8  AP,AiD,All,Bll ,B22,D 

COMMON  FC,  FS,  2PS1  ,  EPS2,  PI,  1  TYPE ,  N,  1TRANS 
COMMQN/F ILT/WC,  W 1 , M2 ,AMP , AO , A 1 , A2 ,  Bl,  B2,  P,  l 
W3=P1*FC/FS 
A2=-1.000 

AP*-castwi«W3)/:as<ui  -M3» 

DO  100  K  =1 *N 

0  =  1 .0D0-81IK  )*AP*B2  IK  J»AP*AP 
AOO=IAO(K)-A1 (K1 * AP >* (1 .ODO-AP I/D 
All  =  tA0(Kl*AP-Al|Kl  1«  Cl.  ODO-AP  1/0 
B11MAP*AP-B1  IK  I  *41  >000«  AP  *AP  H2  >0D0*B2  CK 1  *AP»/D 
B22=  (AP*AP-AP«=31  m  +B2IK  11/D 
AOlK)=AOO 
A1 f K)  =  A1 1 
Bl ( K I s  B 1 1 
B2CK 1  =  822 
100  CONTINUE 
RETURN 
END 


B-8.  SUBROUTINE  FPLOT 


SUBROUTINE  FPLOT 

COMMON  FC,  FS,  EPS1,  EP52,  PI,  1TYPE,  N,  ITRANS 
COMM ON /F ILT/MC,  ¥1 ,M2  ,AMP,AO,AI ,A2,  Bl,  B2,  P,  I 
REA L*8  AMP 

REAL*B  A0C201,  A1I20I,  BH20I,  B2I20),  PI.A2 
DIMENSION  F 12511 ,HE  JMTI251 1, PHI  12511 
DIMENSION  X  l  ABI  2 1 ,YLAB1(2) ,YLAB2(2),HEA0(10),SUBH(2) 
DATA  YLAB2I1 ),YLAB2I2)/8H PHASE  ID  ,BHEG 1  l 

DATA  YLAB1I1 )  >VL  ABI I2)/8HAMPL1TUD,8HE  / 

DATA  XLABI1I.XIABI2 1/ BHF RE QUE NC , 8HV  fF/FSl/ 

VR 1  TEf 6, 800  1 
JJ  = IN* 11/2 
DO  805  I*1,JJ 

8)5  MR  I TE IB, 801 1A0I1 I, All  11,81(11,82(11 

MR1TEI6,802)AMP,A2 

8)0  FORMAT I16X2HA ",32X2HA 1 ,32X2HB1 ,32X2H82 ,/// 1 
801  FORMAT (1 P4 ( 1 XQ32 >25 ) 1 
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B-8.  SUBROUTINE  FPLOT  (Cont'd) 


802  FORMAT I4X6HAMP  =  , I PQ 32 . 25 , 1 OX 5H A?  =  .0PF4.0I 
JJ  = IN* 1 )  /? 

RE  AD  15 , 1 01  i  HEAD  II)»1=1»10) 

13  FORMAT  1 1 0A8  ) 

DO  100  1  *1,250 
A=  1-1 

Fill  -A/249  .<• 

W=2.0*PI *Ft  I  1 
C=CQS(U) 

5  =  5 X  N t  W 1 
PHI  411=0  . 

HE JMTI I ) =0 . 

DO  50  K  =  1 ,  J  J 
CN=A0|K)*A1  |K)*: 

DM*— A 1  IK 1*5 

C  D  =  1 -0  *B 1 1 K ) *C*32 IK !• 42.0*C*C-1 -0 1 
DD  =  B1IK)*S*2.0*32IK  >*C*5 
DD  =  — OD 

EE=CO*CD*DD*DD 
FF=4CN*CO*DN*ODl/EE 
E  E  = I DN*C  D-DC  *CN 1 / EE 
HE JMTI1 ) =HE JUTII  l*FF 
PHI  II) =PH11 1)*EE 
50  CONTINUE 

FF=|1.0D0*A2*CI*HEJHT 111 *5*PH1 I ! ) • A2  *AMP 
EE=ll.0D0*A2*C)*PHI II )-S*HE JNl I I) *A2 
HEJWTtl) =  5QR T IEE  *EE  *F  F  *F  F 1 
PHI  II I  * ATAN2  IEE  , F F  1 

ph  ini  *phii  n*uo.o/p  i 

100  CONTINUE 

CALL  DRAW10I1  ,4,4,20,0,250 ,2.,0.,XLAB, 

1  YLAB1,HEAD, SUBH.F ,HE JWT) 

Fin*o.o 

CALL  DRAN10I1  ,4, 4 ,20 ,0,250  ,2.0,0., XLAB. 
1YLAB2, HEAD, SUBH,F, PHI  I 
RETURN 
END 
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SUBROUTINE  BAN DP 
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SUBROUTINE  BIND? 

CQHHON/F ILT/MC.  U l , M2  ,AMP,A0 ,A1 ,A2 ,  Bit  B2  ,  P,  2 
COMMON  FC,  FS,  EPS1,  EPS 2.  PI,  ITYPE,  N,  ITRANS 
REAL  *8  AOf  20 1  ,  A1I20),  B1I201,  B2I20I,  P>rA?,AMP 
COMPLEX*  16  2t20),P(20),PP»Q,R,S,T  ,U  ,V 
RE  A  L  *8  C  ,D,AL  »E,F 
A2=Q .000 

AL=DC0S(DBLE{Ut*tf2) l/DCO  SIDBIEIW2-R1J) 

JJ-N/2 

IFtN.EQ.DCO  TO  S3 
DO  50  1=1, JJ 
AMP  =AMP * AO (  I  ) 

C=B2fl)~Bltl)**2/4.0D0 
IFIC.LE.O.OOOIGO  TO  40 
PP  =  DCMPLXI-B 1(1} /2 .000 ,DSQRT (C 1 1 

Q=IA0II1*PP*PP*U3< I  MAI <11 I*PP«0CMPLX(A1I1 ) ,0.00011 
l/DCMPLXtO.ODO,2.030*DlMAC(PPl) 

V*DCMPLX<1. 000*3.3001 
S  =  AL*IV«PPl/2  J)30 
R=S*CDSaRT<5*S-PP| 
s=s-cosaRHS*s-ppi 
T*0*  <AL*R-V 1 /  IR“S I 
U=Q*IAL*S— V1/IS-R1 
C  =  DREAL<  R) 

D  =  D  IMAG(R) 

E=DREALIT) 

F  =  D IMAG( T| 

F=-2.0*<C*E*0*F| 

0=C*C*0«0 

C=-2.0*C 

E=2.0*E 

AMP  =AMP*F/D 

AO  (  1  )  =  — F/D 

Alii  l  =  E“F*C/D 

Bill  )sC 

B2(  I  )  =  D 

K=<N*1I/2*I 

C=OREAL(SI 

D=D  IMAG(S) 

E  =  DREALIUI 
F=D IMACIUI 
F=-2.0*IC*E*0*f  I 
0*C*C*0*D 
C  =  -2.0*C 
E  *2 .0*E 
AMP  =AMP«F/0 
A0«RI*-F/0 
A] IK )=E-F*C/0 
B1 IK ) *C 
B2IK)=0 

CO  TO  50  59 
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B-9.  SUBROUTINE  BAN DP  (Cont'd) 


43  C  =  0  SCAT  t  — C  1 

D*-B1II)/2.jD0*C 

C=-B1(1»/2.0D0-C 

EMA0(l)*C*C*(A3tl>«Al(l 1)  *C*A1I  Ill/IOOl 

FMAOI  l)«D*D«(*3ti  MAI  11  ) } *D«A 1 t Ill/ID-C > 

AKP=AMP-E/C-F/D 

AO  I  1 1-E/ C 

All  I 1=-AL*E/C 

B 1 ( 1 )=~AL*( 1 .0D3*C) 

B2(I1*C 
K*IN+1 1/2*1 
A0(K  1  =  F/D 
A1(K)S— AL*F/D 
B1«KJ=-AL*I1. 003*01 
B2 1 K  )*D 
53  CONTINUE 
K*N-2»JJ 

I F IK .EQ.01GG  TO  130 
AMP  =AMP« AO (N)/BilNI 
D-B 1  IN ) 

C*A0INI*(1.CD0-1.3/01 
K  *  I  N*1 1/2 
AO  IK  MC 
A1I K )S-AL*C 
B1(K)*AL •ID-1. 0301 
8  2 1 K 1  *  -D 
130  N=2*N 

RETURN 
END 


B-10.  SUBROUTINE  BANDS 


SUBROUTINE  BANDS 

COMMON/F ILT/MC,  Ml ,M2,AHP,A0,A1»A2,  Bl,  B2.  P,  Z 
COHHON  FC,  FS.  EPSlt  EPS 2»  PI,  HYPE,  N,  ITRANS 
REAL *8  AO (20 )  ,  A1I201  ,  81(201,  62(201,  P I,A2,ANP 
COMPLEX*  16  Z(201,PI201,PP,QtRfS,T  ,U,V 
REALMS  C,D,AL,E,F,X,C1,C2 
A2*0 .0D0 

AL  *DC0S 1 DBLE f  Ml  *M2 1 1/DC05(DBLE (M2-M1 1 1 
JJ*N/2 

X*DTAN(DBLE IM2-M1»I*DTANIP1*FC/FS» 

C1*2.0*AL/(X«1.3D3I 

C2* I 1 .ODO-X 1 / (X* 1 .0  DO  1 
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B-10.  SUBROUTINE  BANDS  (Cont’d) 


IFfN.EQ.JlGC  TO  5 ) 

DO  50 

v*dcmplxh.odo.o.ddoi 

AHP  =  AHP*  AO (  I  ) 

C*B2tl  1-Bdl  )*«2/4-0D0 
1FIC.LE. 0. 000)03  10  AO 
PP=DCHPLXI-B1  I1)/?.0D0»D SORT  ICM 

Q'(AOIJ)  *PP*PP*IA0t  ))  *A1  I!  »l«PP*DC»PLXIAlin  r0 .00011 
1/DCMPl X I  0. ODO, 2. 0D0*D1MAGIPP II 
Q=Q*C2/IV-PP*C2) 

AMP=AHP*2.0*OREAL10) 

S  =  C1«IV-PP)/U¥-C2*PP  1*2 .0 1 
T  =  IC 2*V-PP J/(V-C2*PP) 

R=S*C0SQRTIS*S— T  ) 

S  =  S-C0SQPTIS«=S-f  J 

T=Q*«R*A-Cl/C2***V/C2l/tA-S) 

U=Q*IS*S-tl/C2*S*Y/C21/IS-Rl 

ODREALlRl 

0=0 IMAG ( R ) 

E  =DREAl IT) 
f-DIMAGIT) 

F  =-2.0*IC*F  *0*f  ) 

0=C*C*D*0 
C=-2.0*C 
E  =  2.0*E 
AHP  =AHP*F/D 
A0I 1»*-F/D 
AKI  1=E-F*C/D 
B1 C  1  )- C 
B2I I  )  =  D 
A  =  I N*1  1/2*1 
C  *DREAL IS) 

0=0 1HAGIS1 
E  =DREAL  tU) 

F*0IHAG«U1 
F  =  -2..0*IOE*D*F) 

d=c*c*o«o 

C  =  -2  *0*C 
E*2 .0*E 
AHP  *AHP*F/D 
ADIK)=-F/0 
AllKI=E-F*C/0 
81 IK  )=C 
82IK )=0 
GO  TO  50 


\ 
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B-10.  SUBROUTINE  BANDS  (Cont'd) 


43  C  =D  SORT  ft  — C  I 

D*-Bl«J)/2.0Dd*c 
C  =  -B1I  D/2.0D0-C 

EMAOtl  )*C*C*  IA3C1  MAI  Cl  )l*C«AI  l I ) >/l  C-01 
fHA0m*0*DH*D{I)*AlMH«D«M(I))/ID-C) 
AMP=AMP*E/IC2-C)*F/  1C  2-0  1 
AO  I  l)=E»tC2/fl.3D3«C2)-l  .0/IC2-C)) 

A1 «!)=-£ *C1*(C*C2>/  II .0D0«C2)**2 
Bll  1)  =  C1«(C- 1  .030) /II  -0D0*C 2 ) 

B2( I »  =  €  C  2-C )/fl.0D0«C2) 

K=1*N 

B2IK)-IC  2-D )/tl.0D0«C2) 

B1 IK)  =  C1*(D-1  . 0001/11  .0D0«C2) 

AO  (K  )=F*IC2/ 1  1.300*  C2  1-1  -0/IC2-DN 
A 1  IK)=-F*C1MD«C21/  II  .0D0«C2  1**2 
50  CONTINUE 
K=M-2*JJ 

1FIK.EQ.0)  GO  13  100 
K  =  I  N*1  1/2 
AMP=AMP«AO|K) 

D=B1  IK) 

C*A0CK)/I1.0D0«81IK 1 *C2) *1 1 -ODO-O ) 
B2IK)=IC2*D  )  /  11 . 3  30  *C 2*0  1 
B1  IK)X— C1*ID«1. 3031/11. 000 *C?*D I 
A HP  =AMP«C/B2 (Kl 
A0IK)=C*C2— C/B2IK) 

A 1  (K  ) *— C  1*C -  CM1IH  1/B2IK) 

100  N*2«N 
RETURN 
END 


B- 11.  SUBROUTINE  FILIMP 


SUBROUTINE  FILHP 

COMMON  FC,  FS,  EPS1  ,  EP52,  PI,  I  TYPE ,  N,  1  TRANS 
REAL*8  AMP 

REAl*B  AO  1 20 )  *  AH20),  B1I20),  B2I201,  PI,A2 
COMMON/ VDATA/  VI  N 1 1 00 0 1 , VOUT  11000), NT IMES 
COMHON/F 1LT/NC,  M1,W2,AMP,A0,A1,A2,  Bl,  B2,  P,  Z 
RE A  L*B  YM2.201,  YSUM1,  YSUM2,  YTEM 
YSUM1*0.G 
YSUM2s0 . 0 

jj-(N«n/2 

00  10  K  =  1  »  J  J 
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B-11.  SUBROUTINE  FILIKP  (Cont'd) 


VK(2,K)=A0<K)*VIN(1 > 

VKI1  ,K)  =  A0CK>*tfIN(2HIAl(K)-BllK)»A0(K>)«VlN<l> 
YSUH2sYSUM2+VK(2,K) 

VSUH1=YSUH1*YK41 ,<> 

13  CONTINUE 

VOUT  (1)= Y5UM2 
VQUT 12 ) *  YSUH 1 
DO  50  J)*3,NTIM£S 
Y$UM2*VSUMl 

rsuNi*o.o 

DO  40  K *  1  *  JJ 

YTEM=AO(Kl*VlNtll*A 1 IK 1*V1N IJ-1 1 
1-B1 (K)«YK(1 ,Ki-92IK )*YKI2,K) 

yku,ki«yrii  ,ki 

YK 1 1  ,K  1  *  Y1  EM 
YSUM1*YSUM1*YTE1 
43  CONTINUE 

VQUT  I J)  -Y5UM  1-»YSUN2  *A2  ♦AMP«V!NIJ) 

53  CONTINUE 

RETURN 
END 


B-12.  SUBROUTINE  PLTDAT 


SUBROUTINE  Pt.  TO  AT 

CONMON/VDAT A/VlNt 1 300  1  .YOU! 1 1000 1 , NT  I  ME S 
CONNON  FC.FS 

DIMENSION  HE  ADUO  )«  XL  ABI2) ,VLABI 2|  .TC1000I 
DATA  XLABU>,XLABt2»/8NTIME  iSE ( 8HC0NDS 1  / 

OATA  SUB1.SUB2/SHV1N  .BHVOUT  / 

DATA  YLABI1) , YL AB 1 2 1/BHAMPL 1TUD ,BHE  / 

00  20  1*1»NT 1MES 
23  T( 1 1 *( 1- 1 l/F  S 

RE AO (5,10) ( HEAD* 1), I- 1,101 
13  FORMAT  1 1 OAB  1 

CALL  0RAN10I1  ,4,4 ,20, 2tN  TIME  5, 2.,0  «»XUB,  YLAB ,  HE  AD  ,SUB  1  ,T  ,  V1N 1 
READ (5, 101 (HE ADI  11,1  =  1,101 

CALL  ORAVlOU  ,4, 4 ,20,2  tNT!ME  5,2.,  0.,XL4B»YLAB,  HE  AO,  SUB  2.T,  VOUT’ 

RETURN 

END 
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JOINT  CHIEFS  OF  STAFF 
ATTN  J-3 

WASHINGTON,  IX’  20301 
DIRECTOR 

JOINT  STRATEGIC  TARGET  PLANNING 
STAFF,  JC:: 

ATTN  JSAS 
ATTN  JPST 

ATTN  NR  I -ST  INK')  LIBRARY 
OFFUTT  AFB 
OMAHA,  NB  OHM? 


CHIEF 

LIVERMORE  DIVISION,  FIELD  COMMAND 
DMA 

LAWRENCE  I  I  VERM' (RE  IABORAT'.*RY 
PO  BOX  HUH 
ATTN  P  PRI. 

LIVERMORE,  <  A  #4*-Su 

NATIONAL  '  ‘'MM1 /NIC AT  I*  ’NS  SYSTFM 
OFFICE  OF  THE  MANAGER 
DEPARTMENT  of  f i KEENS E 
ATTN  Nrs-TS,  "MAPLES  I).  BoDSoN 
WASH  I  NOT*  iN,  I*‘  ?U3»j6 
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DIRECTOR 

NATIONAL  SECURITY  AGENCY 
DEPARTMENT  OF  DEFENSE 
ATTN  R-52,  O.  VAN  GUNTEN 
ATTN  S232,  D.  VINCENT 
FT  GEORGE  G.  MEAD,  MD  20755 

UNDER  SECRETARY  OF  DEF  FOR  RSCH  & 
ENGRG 

DEPARTMENT  OF  DEFENSE 
ATTN  G.  BARSE 
ATTN  S&SS  (OS) 

WASHINGTON,  DC  20301 

COMMANDER 

BMD  SYSTEM  COMMAND 
DEPARTMENT  OF  THE  ARMY 
PO  BOX  1500 
ATTN  BMDSC-A OLIB 
HUNTSVILLE,  AL  35807 

COMMANDER 

ERADCOM  TECHNICAL  SUPPORT  ACTIVITY 

DEPARTMENT  OF  THE  ARMY 

ATTN  DRDCO-COM-ME,  G.  GAIJLE 

ATTN  TECHNICAL  LIBRARY  DIV 

ATTN  DELCS-K,  A.  COHEN 

ATTN  DELET-IR,  E.  HUNTER 

FT.  MONMOUTH,  NJ  07703 

COMMANDER 

US  ARMY  ARMOR  CENTER 
ATTN  TECHNICAL  LIBRARY 
FT  KNOX,  KY  4<H21 

COMMANDER 

US  ARMY  COMM -KI.KC  KNGRG  INSTAL 
AGENCY 

ATTN  nCC-PRSO-S 

ATTN  CLX:-CEI)-SF,S 

FT  HUACHUCA,  AZ  86M  1 

COMMANDER 

US  ARMY  COMMUNICATIONS  COMMAND 
ATTN  CX*-OPS-PD 
ATTN  (.V -OPS -OS 
FT  HUACHUCA  AZ  H56 J 3 

COMMANDER 

US  ARMY  COMMUNI"ATI<  »NS  COMMAND 
COMBAT  DEVF..OWKNT  DIVIsD'N 
ATTN  ATS  I  -  ’D-MD 
FT  HUACHtICA,  A/.  HM.I  3 

CHIEF 

US  ARMY  "<  WMUN t"ATl"NS  SYv  AGEN- '  V 
ATTN  CCM-PD-T,  CCM-AD-SV 
FT  MONMOUTH,  N.T  ()77i'3 

PROJBT  OFFICER 
us  A^MY  !  ■;  MM  UN  I  AT  I  "NS  RES 
A.  DEV  "OMMAN1 
ATTN  DR" I'M -A 7”' 

ATTN  f  iRt  ‘  I’M  -Tl  #’>  -BS  1 
FT  M-  )NMf  "TM,  NJ  <i//ni 

D  I  VIS  I  ’  >N  ENG  INFER 

US  ARMY  ENGINEER  DIV,  HUNTS  VI  i.LE 

fu  H"X  iwm,  WEST  STATION 

ATTN  HNI i ED -SR 

ATTN  A.  7.  H"LT 

HUNT "VI I .LE  ,  AL 
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IIS  ARMY  INTEL  THREAT 
ANALYSIS  DETACHMENT 
ROOM  JOI,  BLDG  A 
ARLINGTON  HALL  STATION 
ATTN  RM  7200,  BLDG  A 
ARLINGTON,  VA  222 1J 

COMMANDER 

US  ARMY  INTELLIGENCE  4.  SEC  iW> 
ARLINGTON  HALL  STATION 
4000  ARLINGTON  BLVI) 

Arm  TECHNICAL  LI  HR  ARY 
ATTN  TECH  INFO  FAC 
ARLINGTON,  VA  22212 

•OMMANDER 

IIS  ARMY  MISSILE  COMMAND 

ATTN  DRCFM-PE-EA,  WALLACE  O.  WAGNER 

ATTN  DROPM-PE-EG,  WILLIAM  R.  JOHNSON 

ATTN  DRDMI-TBD 

ATTN  DRDMI-EAA 

REDSTONE  ARSENAL,  A!,  3N80H 

COMMANDER 

US  ARMY  TEST  AND  EVALUATION  COMMAND 
A  H'S  DESTF-FA 

ABERDEEN  PROVING  GROUND,  MD  2  IPOS 
v'i  MMANDER 

US  ARMY  TRAIN  INC,  AND  l  *»  X  'TR  INK 
COMMAND 

ATTN  ATORI-OP-SW 
FT  MONROE,  VA  2  3661 

Ci  MMANDER 

WHITE  SANDS  MISSILE  RANGE 
DEPARTMENT  'K  THE  ARMY 
ATTN  STEWS  -TK -  AN ,  J.  -'K'WA 
WHITE  SANDS  MISSILE  KANoK,  NM  88002 

•  'FED -KR- IN -CHARGE 
CIVIL  ENGINEERING  LABORATORY 
NAVAL  CONSTRUCTION  BATTALD'N  CENTER 
A'TI’N  CODE  I.UHA  (LIBRARY) 

ATTN  CODE  1.08A 

P"RT  MHKNH4K,  1  "A  304  1 

'■  'MMANDER 

NAVAL  AIR  SYSTEMS  COMMAND 
ATTN  AIR-1SOE 
WASHINGTON,  DC  21  3».t> 

-  •’  MMANDER 

NA  VA  f .  KI .  M. '  TR<  'NIC  S  YSTKMS  '  '■  *MM  AN! 1 
ATTN  IMF  11  7-,» IS 
WASHINGTON,  DC  7'Hf.O 

«  1  MMANDER 

NAVAL  I'KAN  SYSTEMS  CENTER 
AITN  CODE  ot-.,  ■.  ELE'P  HER 

’  T  I N  C'-I'E  /24>>,  s.  w.  LICHmAN 
AN  i»1Ei.t«,  -  A  ■*.’  1  H2 

COMMANDING  OI'FICER 
NAVAL  oRDNAf*  l  TATD'N 
ATTN  S  I'ANDARD  I  .’.AT  D  f.  '  IV 
INI 'I  AN  HEAD,  MD  70M 

sobeRINTENDKNT  (CODE  1424) 

NAVAL  POSTGRADUATE  SCHooI, 

A  ITN  ■  •-  'I»E  1424 
M-'NTEREY,  'A  •»  *,  *4,1 


D I R  ECTOR 

NAVAL  RESEARCH  LABORATORY 

ATTN  CODE  4104,  WANSAL  L.  BRANCATo 

ATTN  CODE  2t>27,  U*R  IS  R.  FOI.EN 

ATTN  CODE  St, 21,  RICHARD  L.  STATLER 

ATTN  CODE  *>*,24 

WAS H I NGTON ,  fX’  20  i  7  6 

COMMANDER 

NAVAL  SHIP  ENGINEERING  CENTER 
DEPARTMENT  OF  THE  NAVY 
ATTN  CODE  6174P2,  EIHVARI'  F.  DUFFY 
WASHINGTON,  IX'  20  362 

COMMANDER 

NAVAL  SURFACE  WEAPONS  CENTER 
ATTN  CODE  WA51RH,  RM  130-108 
ATTN  CODE  F32,  EDWIN  K.  RATHBURN 
WHITE  OAK,  SILVER  SPRING,  MD  20'Mi‘ 

COMMANDER 

NAVAL  SURFACE  WEAPONS  CENTER 
DAHLGREN  lABORATORy 
ATTN  CODE  DF-56 
DAHLGREN,  VA  7244s 

v 'oMMANDKR 

NAVAL  WEAPONS  CENTER 
ATTN  CODE  S3 3,  TFX’H  LIB 
CHINA  LAKE,  CA  U3SSS 

COMMANDING  OFFICER 

NAVAL  WEAPON: S'  EVALUAT'D >N  FACILITY 
KIRTIAND  AIR  FORCE  RASE 
ATTN  CODE  AT-b 
ALBU^OER-JOK,  NM  8711/ 

OFFICE  OK  NAVAL  RESEARCH 

ATTN  CODE  427 

ARI  1NGT>'N,  VA  2721.* 

DIME*  "IV'R 

STRATEGIC  SYSTEMS  PROJECT  OFFICE 
NAVY  DEPARTMENT 

ATTN  NS p- 770!  ,  JOHN  W.  ITTSKN8KKGFK 

ATTN  NSP-2M7,  RICHARD  L.  COLEMAN 

ATTN  NSp-4  3,  TECH  I  IB 

ATTN  NSp  - 7  M  <4 

ATTN  NSP-7  30,  |i.  u  'Ll 

WASH  INGT*  >N,  DC  2«M7b 

•  MMANDKh 

AERONA  TD’AL  SYSTEMS  DIVISION,  AKSi 
A  rm  ASD-YH-FX 
ATTN  KNK  TV 
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PATRICK  AFB ,  EL  3  7  'L’‘> 
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Arm  ntn 

A  TTN  NT 

ATTN  EL,  -  ARI.  E.  BAUM 
A  ITN  EL  XT 
ATTN  SIT. 

Arm  ca 

ATTN  El  A,  I  .  'AST  ILL  > 

A  TTN  ELI 

ATTN  ELT,  W.  PAGE 
A'fTN  NXs 

KIRTLANI  AFP ,  KM  M  ’  I  1  ' 


DIRECTOR 

AIR  UNIVERSITY  LIBRARY 
DEPARTMENT  OF  THE  AIR  FORCE 
AITN  AUL-LRK-70-2S0 
MAXWELL  AFB,  AL  3b1  1  2 

HKADOCARTKkS 

ELECTRONIC  SYSTEMS  DI VIS  ION/ YS KA 
DEPARTMENT  OF  THE  AIR  FORCE 
ATTN  YS  EA 

HANSCOM  AFB,  MA  01731 
C'TMMANDER 

FOREIGN  TFX'HNODIGY  DIVISION,  AFSC 
A'lTN  NICD  LIBRARY 
ATTN  EXIT,  R.  L.  KALIAKD 
WRIGHT -PATTERSON  AKR,  »H  4^4  33 
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1  >Gl > FM  AIL  MM  EDGE 
DEPARTMENT  OF  THE  A.R  KORCR; 

ATTN  C-A1L  MMKTH ,  P.  W.  HER  THE! 
ATTN  MM ED",  LHC  KIDMAN 
ATTN  MAJ  R.  BIACKBUKN 
HILL  AFP,  iT  R44"t. 

Ci.  MMANDER 

Rv-WK  AIR  DEVELOPMENT  CENTER,  AFSC 
A  TIM  TSLU 

v,R  IFF  ISS  AFB,  NY  13441 
Ci  MMANDKH 

SA<  'HAMHh'T* '  AIK  Is  JGJSTIt  *S  i 'ENTER 
DEPARTMENT  OF  THE  AIR  FORCE 
ATTN  MMCRS,  H.  A.  PKIMASTK' 

AT’IM  MM  IRA,  .1.  W.  DEMES 
ATTN  MMSREM,  F.  K.  SPEAK 
MCCLELLAN  AKR,  CA  '*Sb‘0 

SAM So  IN 

AIR  FORCE  SYsTEMS  C  'MM AND 
po  BOX  '*7 •*«,«' 

WORLDWAY  POSTAL  CENTER 
(  INTELLIGENCE) 

A 1TN  INI' 

LOS  ANGELES,  CA 

SAM  Si '  MN 

AIR  FORCE  SYSTEMS  •  MMAND 
(  M  l  N  ITEM  AN  ) 

AITN  MNNH,  MAS  M.  BAR  AN 

ATTN  MNNH,  .'AIT  R.  LAWREN'  I 

N.'RT'N  AFB,  CA  »7  4"  * 

sAMS  •  YA 

A  I  F'  'Rt ' F  S  Y'-T KM:  '•  MM  AND 
D *  B<'X  •»7,«iH' 

W  iRLlWAY  POSTAL  ENTER 

Arm  ya  (t  ■ 

’..os  ANGKLKS,  CA  Hi’i'H" 

STR  A  I’F'  -  I1'  MR  .  WMAND  XI  I  s 
ATTN  NR  I  -  ,TINK  ’  !  I  PR  ARY 
ATTN  DEL 

ATTN  <iARN!"l  I.  MAT'S  K I 

AITN  XPK  • ,  MAS  BRIAN  G.  •  H  PHAN 

■  lKi-'.'TT  AC"  ,  NR  •  1  1  3 

!  FPAR  IMI  N'l  •  INERCY 
A 1 1"  •>>(•>■. ,'j'f  TER  At  I  'N.  'If  I  >' 

D  ■  B">  ‘->4".' 

AITN  VFU  I  I  BEAR  V 

Aim  HFRATIonM  ;>  A '  I  TV  D1V 

A!.P--.,,!ER  .  IT  ,  NM  »*'u  i «, 
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UNIVERSITY  OF  CALIFORNIA 

LAWRENCE  LIVERMORE  LABORATORY 

PC  BOX  808 

ATTN  THOH  INFO  DEPT 

AITN  L-Ob,  T.  :k>NICH 

ATTN  L-S4*>,  D.  MEEKER 

AITN  1.  -  I  Sb,  K.  MILLER 

AITN  1.-10,  H.  KROGER 

ATTN  H.  S.  OAHAYAN 

LIVERMORE,  CA  94SS0 

LOS  ALAMOS  SCIENTIFIC  LABORATORY 

PO  Ht>X  IbbJ 

ATTO  BRUCE  W.  NOEL 

ATTN  CLARENCE  BENTON 

LOS  ALAMOS,  W  87S4S 

SAND I A  LABORATORIES 
LX'  HO.X  SBOO 
ATTN  C.  N.  V  ITT  I 'POE 
ATTN  R.  L.  PARKER 
ATTN  ELMER  K.  HARTMAN 
albucukrcue,  nm  b7hs 

CENTRAL  INTELL.  IGENCK  AGENCY 
Arm  RP  'SI,  RM  bC.48,  HC,  HI. Pi  I 
(OS  I  NKP  'NWB 1 
WASHINGTON,  DO  20b0b 

A ’PM  IN  IS  THAT'  >R 

DEFENSE  ELECTRIC  POWER  APMIN 
DEPARTMENT  OF  THE  INTERIOR 
INTERIOR  SOUTH  BLtK,! ,  112 

ATTN  L.  o*NEU,I. 

WASHINGTON,  lx-  J0J4O 

i  'KPAkTMKNT  OF  TRANSPORTATION 
FEDERAL  AVIATION  ADMINISTRATION 
HEAP-,  FARTERS  SEV  PIV,  ASK- .10*' 

Mi ’* 1  INDEPENDENCE  AVENUE,  SW 
ATTN  NEC  HIV  ASF  -  300 
WASHINGTON,  if  .’0S'»1 

AEROSPACE  CORPORATION 
P"  lii  >x 

Arm  C,  B.  PF.ARI.STON 
ATTN  IHVING  M.  GARFPNKFL 
ATTN  .LILIAN  RKINHFIMER 
ATTN  LIBRARY 
ATTN  CHARLES  GRKFNHOW 
(.OS  ANGELES,  GA  *1000* 1 

A  LAB  I  AN  ASSiX  •  IATES 
."'P  NORTH  NASH  STREET 
ATTN  LIBRARY 
EL  SECOND* >,  oa  *in.,4S 

AVI  ‘  RESEARCH  k  SYSTEMS  GRoMP 
.*01  DWELL  STREET 
WttMINGToN,  MA  olHM/ 

BATTEL;  F  MWORIAL  INST  IT'  :TI 
Si-S  KIN,,  AVKNI'F 
ATTN  ROBERT  H.  BI.A2.KK 
ATTN  El -GENE  R.  LEACH 
i  *■  *•  l MBPS ,  on  4  l.’OI 

H«  ’*  )R  L**  >HAT  I ■  *N 

■*1  .  .MNEs  BRANCH  DRIVE 

MTN  •  »Rpi  iRATE  LIBRARY 
MO!. I  AN,  VA  2/»  V 


BUM  CORPORATION 
VO  BOX  9274 

ALBIVPERUDE  INTERNATIONA!. 

ATTN  LIB 

AJ.BUCl'EROPE,  W  87119 

BENnIX  CORPORATION 

RESEARi ‘ H  I ABORATX  'R  I ES  f  I VI S  I « >N 

BENI > IX  CENTER 

ATTN  MAX  FRANK 

SOUTHFIELD,  MI  4807  S 

BENPIX  CORPORATliW 
NAVIGATION  AND  CONTROL  GR-*:;L 
ATTN  PEPT  MP1 
TKTERRORO,  NO  07*08 

BOEING  COMPANY 
PO  BOX  3707 

AITN  Hi VIARP  W.  WIPKLEIN 
ATTN  P.  K.  ISBFLI. 

ATTN  PA  VIP  KFMI.E 
ATTN  h.  C.  HANRAHAN 
ATTN  KENT  TECH  LIB 
SEATTLE,  WA  98124 


HRiWN  ENG  INKER  I NG  COMPANY,  INC 
CIMMINGS  RESEARCH  PARK 
ATTN  FRF.P  L  EON  ARP 
HUNTSVILLE,  AI.  3S8H7 


BURROUGHS  CORPORATION 
EKPERA1.  AN!'  SPECIAL  SYSTEMS  GROI'P 
CENTRAL  AVF  ANI'  HOUTE  2S2 
D'  BOX  b*  1  7 

ATTN  ANGELO  J.  MAURIKLIG 

PAoi.i,  pa  t  noi 


CALS PAN  \  'R  IN  'RAT  Ii  'N 
LX'  Box  4 HO 
ATTN  TECH  LIBRARY 
BUFFALO,  NY  M.V1* 


CHARLES  STARK  PR A PER  lAHORATi'H Y , 
INI'. 

S‘>S  TECHNOLOGY  SvHIARF 
ATTN  KENNETH  FKRT1G 
ATTN  Tio  MS  14 
CAMRR  I  t'GK,  MA  O.M  .V* 


CINCINNATI  ELECTRON  B  *b  Ci  'K  pi  'RATION 
.'Bio  t  ;i  ENf  -ALE  -  Mil  -K.  'V,  R*  'AI 
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ATTN  M  V  IN  >•  H  I  FF 
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CONTROL  PATA  CORPORATION 
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ATTN  JACK  MEEHAN 
MINNEAPOLIS,  MN  bS440 

CUTLER -HAMMER,  INC. 

All.  DIVISION 
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PEEP  PARK,  NY  1  1  72'1 

PIKEWDOP  CORPORATION 
1  M  l  UNIVERSITY  BI.VP,  NF 
ATTN  TECH  LIB 
ATTN  L.  WAYNE  DAVIS 
ALBU'jUER  >,'CF,  NM  87102 


DlKEWiOD  CORPORATION 

271 1,  i V KAN  &  PARK  BI.VP,  SUITE 

ATTN  K.  LEE 

SANTA  M'NICA,  CA 


E -SYSTEMS,  INl 
GHEKNV 1 1,1  E  PI  VIS  I  'N 
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ATTN  JOi.FTA  MU 'RE 
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